A converse approach to the calculation of NMR shielding tensors 
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We introduce an alternative approach to the first-principles calculation of NMR shielding tensors. 
These are obtained from the derivative of the orbital magnetization with respect to the application 
of a microscopic, localized magnetic dipole. The approach is simple, general, and can be applied to 
either isolated or periodic systems. Calculated results for simple hydrocarbons, crystalline diamond, 
and liquid water show very good agreement with established methods and experimental results. 
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Nuclear magnetic resonance (NMR) measures the tran- 
sition frequencies for the reorientation of nuclear mag- 
netic moments in an applied magnetic field. Since the 
local magnetic field differs from the external one as a 
result of electronic screening, NMR spectroscopy [1] has 
been recognized since 1938 [2\ to be a powerful experi- 
mental probe of local chemical environments, including 
structural and functional information on molecules, liq- 
uids, and increasingly, on solid-state systems. 

First-principles calculations of NMR spectra were first 
developed in the quantum chemistry community Q and 
applied to molecules and clusters, but applications to ex- 
tended crystalline systems were hindered by the difficulty 
of including macroscopic magnetic fields, which require 
a non-periodic vector potential that is not compatible 
with Bloch symmetry. In 1996 Mauri et al. developed a 
linear-response approach for calculating NMR shieldings 
in periodic crystals based on the long- wavelength limit of 
a periodic modulation of the applied magnetic field [J]. 
In principle, this long-wavelength approach could be im- 
plemented using standard ground-state calculations, but 
only at the cost of introducing prohibitively large super- 
cells. In 2001 Sebastiani and Parrinello used a local- 
ized Wannier representation to derive an alternative 
linear-response approach based on the application of an 
infinitesimal uniform magnetic field Q. More recently, 
attention has focused on the development of these ap- 
proaches in the context of pseudopotentials 0, @|, re- 
lying on ideas developed within the projector-augmented 
wave method (PAW) leading to a growing use of these 
methods in combination with modern plane-wave pseu- 
dopotential codes [13, EH- Despite these advances, ex- 
isting methods for computing NMR shifts in crystalline 
systems remain complex, in that they require a linear- 
response implementation with significant extra coding. 

In this Letter, we reformulate the problem of com- 
puting NMR shielding tensors so that the need for a 
linear-response framework is circumvented. For clarity, 
the previous formulations shall be referred to as direct ap- 
proaches, in that a magnetic field is applied and the local 



field at the nucleus is computed. Our alternative, con- 
verse approach obtains the NMR shifts instead from the 
macroscopic magnetization induced by magnetic point 
dipoles placed at the nuclear sites of interest. This ap- 
proach is made possible in periodic systems by the recent 
developments that have led to the Berry-phase modern 



theory of magnetization [Ij, |13|, [ij, [l5|. We demon- 
strate the method by a first application to molecules, 
crystals, and liquids, finding excellent agreement both 
with experiment and with calculations using the direct 
approach. Our new method is simple and general, and 
provides a straightforward alternative avenue to the com- 
putation of NMR shifts, which is suited for large-scale 
simulations and situations where a linear-response for- 
mulation is cumbersome or unfeasible. 

Let us start by considering a sample to which a con- 
stant external mag netic field B"''* is applied. The field 
induces a current that, in turn, induces a magnetic field 
B'"'*(r) such that the total magnetic field is B(r) = 
goxt ^_ Bi"d(i.) In NMR experiments the apphed fields 
are small compared to the typical electronic scales; the 
absolute chemical shielding tensor ct is then defined via 
the linear relationship 
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The index s indicates that the corresponding quantity 
is to be taken at position r^, i.e., the site of nucleus s. 
NMR experiments usually report the isotropic shielding 
as = ^Tr[CTs] via a chemical shift that is defined by 
convention as = Cicf — CTs, where Crof is the isotropic 
shielding of a reference compound. 

As mentioned above, direct approaches 0, S, S S] cal- 
culate the chemical shielding from the current response 
of the system to an external magnetic field, applied us- 
ing perturbation theory and taking the long-wavelength 
limit. The approach we propose is fundamentally dif- 
ferent: instead of determining the current response to a 
magnetic field, we derive chemical shifts from the orbital 
magnetization induced by a magnetic dipole. This can 
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be shown using a thermodynamic relationship between 
mixed partial derivatives [l6|, as follows: Using Bg^a = 
-Ba'' + -Bi"^> Eq. © becomes d^p-a,^^^ = dB,^^/dBf\ 
For the moment, we assume that B°'^' can be replaced by 
the total macroscopic -B-field in the denominator of this 
equation, thus neglecting the macroscopic induced field 
(this restriction, appropriate for normal components in 
a slab geometry, will be relaxed shortly). The numera- 
tor may be written as Bs^a — -^dE/drris^a, where E can 
be interpreted either as the energy of a virtual magnetic 
dipole nis at one nuclear center in the field B for a fi- 
nite system, or as the energy per cell of a periodic lattice 
of such dipoles; we adopt the latter view. Then, writing 
the macroscopic magnetization as Mfs = — dE/dBp 
(where is the cell volume), we obtain 
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Thus, (Ts accounts for the shielding contribution to the 
macroscopic magnetization induced by a magnetic point 
dipole nis sitting at nucleus and all of its periodic repli- 
cas. In other words, instead of applying a constant (or 
long-wavelength) field B'"'* to an infinite periodic sys- 
tem and calculating the induced field at all equivalent 
nuclei s, we apply an infinite array of magnetic dipoles 
to all equivalent sites_s and calculate the change in or- 
bital magnetization 
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Since the pertur- 
bation is now periodic, it can easily be computed us- 
ing finite differences of ground-state calculations. This 
is our principal result. Note that M = nis/ri -|- M'"'^, 
where the first term is present merely because we have 
included magnetic dipoles by hand. It follows that the 
shielding is related to the true induced magnetization via 

It is useful to pause here and consider the analogy with 
the Born [l3| effective charge tensor Z* which may be 
regarded as (i) the component of the force in direction 
a on site by a unit macroscopic electric field £ in direc- 
tion /3 (at zero nuclear displacement), or, alternatively, 
as (ii) the /3-component of the macroscopic electric po- 
larization P linearly induced by a unit displacement of 
nucleus s and its periodic replicas in direction a, in a 
vanishing macroscopic electric field. Since the force on 
nucleus s is given by Fg.a ~ —dE/drg a, (i) and (ii) are 
related by 
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in close analogy with Eq. Note that, in order to 

comply with the Born definition, one must choose the 
lattice-periodical solution of Poisson's equation, corre- 
sponding to vanishing macroscopic electric field. Other 
choices are possible and lead to other kinds of effective 
charges, which are related to each other via the dielec- 
tric constant [l^. By comparing Eq. ^ to Eq. ^ we 



notice that the genuine analogue to is 1 — (and 
not'as), as indeed the names "effective" vs. "shielding" 
imply. As a side note, it should be pointed out that this 
analogy between the electric and magnetic cases would 
be formally more direct if the nucleus carried a magnetic 
monopole charge. 

As in the electrical case [l3|, the choice of magnetic 
boundary conditions implies a choice for the shape of the 
macroscopic finite sample. Following Ref. [l^, shape ef- 
fects can be embedded in the depolarization coefficients 
ria (with J2a ~ 1)' whose special cases are the sphere 
{nx~ny^nz = l/3), the cylinder along z (ux—ny— 1/2, 



-0), and the slab normal to z (n^ 



=0, n^=l). The 



main relationship for the macroscopic fields in Gaussian 
units may be written as Ba — B'j^^ -I- 47r(l — Ua) Ma- 
lt can be seen that for the slab geometry the normal 
component of B coincides with the one of B°'^*. Hence 
our computed crg ^z are suitable for direct comparison 
with measurements of the normal component performed 
on a slab-shaped sample. Assuming non-magnetic media 
with small, isotropic susceptibility it can be shown 
that the shielding for a general shape is related to our 
calculated one by o'f^^p° — '^s,af3 — Sap4:Trx{l — rip). For 
the special case of a spherical sample we have cr^^^°p° ~ 

0-s,a/3 - (87r/3)x<5a/3- 

In order to calculate the shielding tensor of nucleus 
s using Eq. ([2|), it is necessary to calculate the induced 
orbital magnetization due to the presence of an array of 
point magnetic dipoles nis at all equivalent sites . The 
vector potential of a single dipole in Gaussian units is 
given by [lO] 
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For an array of magnetic dipoles A(r) — As(r — R), 
where R is a lattice vector. Since A is periodic, the 
average of its magnetic field V x A over the unit cell 
vanishes; thus, the eigenstates of the Hamiltonian re- 
main Bloch-representable. In the Fourier representation 
A(r) = EG^oA(G)e'G- with 

A(G)^-^^e--^ (5) 

where the reciprocal lattice vector G = may be ex- 
cluded from the sum with no loss of generality. Note that 
we have implicitly chosen the transverse gauge V • A — 0, 
which is apparent from G • (nis x G) = 0. The periodic 
vector potential A(r) can now be included in the Hamil- 
tonian with the usual substitution for the momentum 
operator p ^ p — -A. As a result, the kinetic energy 
operator becomes 
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where nie is the electronic mass and c is the speed of 
light. Due to our choice of gauge, p and A commute. 
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We can now calculate the shielding according to Eq. ([2]) 
by solving for the ground state with the additional terms 
of Eq. ^ included in the Hamiltonian [21[, and then 
calculating the resulting change in orbital magnetization. 

The converse method can be implemented directly 
in any all-electron electronic-structure code. However, 
many popular density-functional theory codes use pseu- 
dopotentials to increase computational efficiency. In or- 
der to calculate NMR shifts in the presence of pseudopo- 
tentials, a PAW reconstruction needs to be performed, 
as shown by Pickard and Mauri We have devel- 

oped this reconstruction methodology for the converse 
method; the rather involved mathematical formalism will 
be presented elsewhere 22]. We implemented our con- 
verse approach, including this reconstruction and the cal- 
culation of orbital magnetization 12, 13, 14, 3], into the 
PWSCF package of the Quantum-ESPRESSO distri- 
bution [23l . We use the PBE exchange-correlation func- 
tional [2 4 and TrouUier-Martins pseudopotentials [2^, 
with convergence of the NMR shifts for a kinetic-energy 
cutoff of 80 Ryd. The dipole perturbation \ms\ used is 
1/LtB, although we find identical results for any value in 
the broad range 10~^ to lO'^/ie- 

As an initial test, we applied our converse approach 
to the calculation of hydrogen NMR chemical shieldings 
for small molecules using a supercell geometry (in such 
cases the full orbital-magnetization theory is not needed, 
as the magnetization can be obtained just by integrating 
the orbital currents over the molecular region, but the 
results presented here were in fact obtained using the 
Berry-phase modern theory of magnetization) . For pur- 
poses of comparison, we also calculated these shieldings 
with the direct method, also implemented by some of us 
in PWsCF, largely eliminating discrepancies due to any 
technicality. The results for these two approaches are 
shown in Table H] together with the experimental values. 
It is immediately obvious that the direct and converse 
methods give almost identical results, validating our ap- 
proach, and also very good agreement (we suspect that 
the slight deviations between the two calculations are 
due to the long-wavelength approximation of the direct 
method) . 

Next, we applied our method to crystalline diamond. 
For our calculations we used an 8-atom cubic cell with a 
lattice constant of 3.498 A. The NMR shielding converged 
to within 0.1 ppm for a k-point mesh of 8 x 8 x 8. For the 
^^C shielding we find 131.20 ppm, in perfect agreement 
with the direct method. To estimate the effect of the spu- 
rious interactions of the localized dipole with its images 
in neighboring supercells, we repeated the calculation for 
a 64-atom cubic cell, finding an almost identical shield- 
ing of 131.24 ppm. The fast convergence with respect to 
supercell size is due to the fast decay (l/r^) of the vector 
potential in real space. 

Finally, we applied the converse approach to com- 
pute the hydrogen chemical shifts in a supercell simu- 



TABLE I: Hydrogen NMR chemical shielding a, in ppm, for 
several diflterent molecules. Structural parameters were taken 
from footnote 22 of Ref. [i]. 





experiment 


direct 


converse 


H2 


26.26" 


26.2 


26.2 


HF 


28.51" 


28.4 


28.5 


CH4 


30.61 " 


30.8 


31.0 


C2H2 


29.26' 


28.8 


28.9 


C2H4 


25.43 


24.7 


24.8 


C2H6 


29.86'' 


30.2 


30.4 



" Reference 
'' Reference 



lation of liquid water. Our supercell contained 64 water 
molecules, twice the size of the largest supercell used in 
previous NMR calculations on liquid water using the di- 
rect method 2^ 2§|. We obtained the atomic trajecto- 
ries from a molecular-dynamics simulation using TIP4P 
[30t potentials under standard conditions. For five snap- 
shots separated by 200 ps, we took the atomic positions 
and thermalized the hydrogen atoms alone for 2 ps using 
ab-initio Car-Parrinello molecular dynamics. This pro- 
cedure is aimed at obtaining a more realistic description 
of the detailed structure of the water molecules while 
retaining the accuracy of the oxygen-oxygen pair corre- 
lation function. We then used these positions to calcu- 
late hydrogen shifts with the converse method. We cal- 
culated the shift of liquid water with reference to the 
gas-phase shift, i.e., iJuquid = CTgas - criiquid, thus re- 
porting the experimental measurable change for the gas- 
liquid transition. For the gas-phase shift the converse 
method gives 31.0 ppm. For the susceptibility correc- 
tion of periodic water we used the experimental value 
for Xwator = —7.2 X 10^^ emu under standard conditions 
[sH . Our distribution for the hydrogen shifts, shown in 
Fig. [U can be directly compared to results obtained us- 
ing the direct method reported in Fig. 4 of Ref. [2^ and 
in Ref. \2^. We find an average shift of 5.94 ppm from 
our distribution, in excellent agreement with 5.83 ppm 
and 5.15 ppm from Refs. 28| and [2§|, respectively. Fur- 
thermore, the spread of our distribution as measured by 
the standard deviation is 2.4 ppm, again in precise agree- 
ment with the value of 2.4 ppm obtained from the direct 
method 

At first sight it might appear that the converse method 
is computationally more demanding than the direct 
method, since we need to perform 3iV calculations to 
obtain the shielding tensor for N atoms. Often, though, 
only a few selected shifts are needed, and even in the 
worst-case scenario (such as the water calculation shown 
before) it should be stressed that re-minimizing the elec- 
tronic wavefunctions in the presence of the perturbation 
is very fast, usually requiring a single self-consistent it- 
eration. This enables the calculation of NMR shielding 
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FIG. 1: Distribution of hydrogen NMR shifts in liquid water 
relative to the gas-phase shift. The distribution was obtained 
from five snapshots of the 64-molecule system (640 hydrogen 
atoms). The dashed line is a polynomial fit and it serves as a 
guide to the eye. The vertical arrow and horizontal line indi- 
cate the position of the average and the range of the standard 
deviation, respectively. 

tensors for systems with several hundred atoms even on 
small computer clusters. However, the main advantage of 
the converse method is the simplicity of its implementa- 
tion, in that it works via finite differences of ground-state 
calculations and does not require a linear-response im- 
plementation. This is likely to be a significant advantage 
for future applications in conjunction with more com- 
plex forms of exchange-correlation functionals such as 
DFT-I-U, exact exchange, hybrid functionals, or beyond- 
DFT correlated-electrons methods. 

In conclusion, we have derived an alternative first- 
principles method for calculating NMR chemical shield- 
ing tensors in solids. The new approach is considerably 
simpler than existing techniques, avoiding altogether the 
difficulties related to the choice of a gauge origin and 
the need for a linear-response implementation. We have 
demonstrated the correctness and viability of our ap- 
proach by calculating chemical shieldings in isolated and 
periodic systems: simple molecules, crystalline diamond, 
and liquid water, finding excellent agreement with pre- 
vious theoretical and experimental results. Applications 
to more complex systems, including small proteins and 
systems with heavy elements, are currently in progress. 
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